Surface pattern formation and scaling described by conserved lattice gases 
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We extend our 2 + 1 dimensional discrete growth model (PRE 79, 021125 (2009)) with con- 
served, local exchange dynamics of octahedra, describing surface diffusion. A roughening process 
was realized by uphill diffusion and curvature dependence. By mapping the slopes onto particles 
two-dimensional, nonequilibrium binary lattice model emerge, in which the (smoothing/roughening) 
surface diffusion can be described by attracting or repelling motion of oriented dimers. The binary 
representation allows simulations on very large size and time scales. We provide numerical evidence 
for Mullins-Herring or molecular beam epitaxy class scaling of the surface width. The competition 
of inverse Mullins-Herring diffusion with a smoothing deposition, which corresponds to a Kardar- 
Parisi-Zhang (KPZ) process generates different patterns: dots or ripples. We analyze numerically 
the scaling and wavelength growth behavior in these models. In particular we confirm by large size 
simulations that the KPZ type of scaling is stable against the addition of this surface diffusion, 
hence this is the asymptotic behavior of the Kuramoto-Sivashinsky equation as conjectured by field 
theory in two dimensions, but has been debated numerically. If very strong, normal surface diffusion 
is added to a KPZ process we observe smooth surfaces with logarithmic growth, which can describe 
the mean-field behavior of the strong-coupling KPZ class. We show that ripple coarsening occurs if 
parallel surface currents are present, otherwise logarithmic behavior emerges. 

PACS numbers: 05.70.Ln, 05.70.Np, 82.20.Wt 



I. INTRODUCTION 



In nanotechnologies large areas of nanopatterns are 
needed, which can be fabricated today only by expen- 
sive techniques, e.g. electron beam lithography or di- 
rect writing with electron and ion beams. Besides the 
conventional 'top-down' technologies, which use masks, 
photoresits, ... etc. to create structures on the surfaces, 
nowadays 'bottom-up' approaches are getting close to 
achieve the same results more efficiently. In that case 
the self-assembly of patterns of large areas is facilitated 
in a cost effective way [l|. This has led reopening of 
the research for fundamental theoretical understanding 
of the ion-beam induced surface patterning and scaling 
Q , which was flourishing at the end of the previous cen- 
tury [H, |H . Although the basic universality classes and 
important models have been explored, many notoriously 
difficult fundamental questions have been unanswered. 
Perturbative renormalization group methods and analyt- 
ical tools have limited applicability and precise numeri- 
cal simulations, approaching asymptotic scaling regimes 
were feasible in one dimension mainly. 

One of the most fundamental problem of kinetic rough- 
ening can be characterized by the Kardar-Parisi-Zhang 
(KPZ) equation [5|. The KPZ has been found to describe 
other important physical phenomena such as randomly 
stirred fluid @, dissipative transport 7, 8], directed poly- 
mers Q and the magnetic flux lines in superconductors 
10 . Therefore, we started our studies, by setting up 
the simplest possible microscopic model exhibiting this 
behavior [llEl. 



The KPZ is a non-linear stochastic differential equa- 
tion, describing the dynamics of growth processes in the 
thermodynamic limit specified by the height function 
h(x,t) 

d t h(x,t) = u + crV 2 /i(x,t) + A 2 (V/i(x,t)) 2 +r/(x,t) . (1) 

Here v and A2 are the amplitudes of the mean and local 
growth velocities respectively, a is a smoothing surface 
tension coefficient and 77 roughens the surface by a zero- 
average, Gaussian noise field exhibiting the variance 



( V (x,t)r ] (x',t , ))=2DS d (x-x , )(t-t') 



(2) 



We denote the spatial dimensions of the surface by d 
and the noise amplitude by D. The pure KPZ equa- 
tion is exactly solvable in 1 + Id [9(, but in higher di- 
mensions only approximate solutions are available (see 
U). In d > 1 spatial dimensions due to the competition 
between the roughening and smoothing, models charac- 
terized by {1} , exhibit a roughening phase transition be- 
tween a weak-coupling regime (A2 < AjjS), governed by the 
A2 = Edwards- Wilkinson (EW) fixed point [13j . and a 
strong coupling phase. The strong coupling fixed point 
is inaccessible by perturbative renormalization method. 
Therefore, the KPZ phase space has been the subject of 
controversies and the value of the upper critical dimen- 
sion is an active field of studies for a long time. 

Mapping of surface growth onto reaction-diffusion sys- 
tem allows effective numerical simulations and better un- 
derstanding of basic universality classes [p34l6l ]. The 
principal aim of this paper is to show that some of most 
fundamental growth processes can be well described by 
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the simplest, restricted solid on solid (RSOS) model with 
Ah = ±1. This strong condition enables a mapping onto 
binary lattice gases and facilitates to create fast algo- 
rithms. We will discuss models, which follow universal 
scaling laws and exhibit pattern formation. Although 
the understanding of coarsening phenomena [171 ] of sur- 
face patterns has been developing rapidly by continuum 
approaches [l8| and the agreement with the ion beam 
induced nanopattern experiments is improving (ljjj . the 
identification of the various coarsening scenarios is still a 
theoretical challenge [20j . Our approach is based on the 
KPZ models we introduced very recently [111, EH, hence 
for the sake of completeness we review them now. 

In one dimension a discrete RSOS realization of the 
KPZ growth is equivalent [HI, HI! to the asymmetric sim- 
ple exclusion process (ASEP) of particles [23[, while we 
have shown that this 'roof-top model' can be generalized 
to higher dimensions [ill, El] • This mapping is interest- 
ing not conceptually only, linking nonequilibrium surface 
growth with the dynamics of driven lattice gases [13, HE] , 
but provides an efficient numerical simulation tool for 
investigating debated and unresolved problems. 

The surface built up from the octahedra can be rep- 
resented by the edges meeting in the up/down middle 
vertexes (see Fig. [1J. The up edges, in \ — x or 




FIG. 1: (Color online) Mapping of the 2 + 1 dimensional sur- 
face diffusion onto a 2d particle model (bullets). Surface at- 
tachment (with probability p) and detachment (with proba- 
bility q) corresponds to Kawasaki exchanges of particles, or 
to anisotropic migration of dimers in the bisectrix direction 
of the x and y axes. The crossing points of dashed lines 
show the base sub-lattice to be updated. Thick solid line on 
the surface shows the y cross-section, reminding us to the 
one-dimensional roof-top model. When the shown desorp- 
tion/absorption steps are executed simultaneously, they real- 
ize a surface diffusion step of size s — 3 along the y axis. This 
corresponds to a pair of dimer movement, a repulsion in case 
of this smoothing reaction. 



X = V directions at the lattice site i,j are represented 
by a x (i,j) = +1, while the down ones by <r x (i,j) = —1 
slopes. Therefore, we approximate surfaces using RSOS 
model with Ah = ±1. 
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In this paper we show that one can describe various, 
more complex surface process without the need of having 
larger Ah height differences. These are built up from 
the basic octahedron deposition/removal processes 
Let us remind the reader that a single site deposition 
flips four edges, which means two , +V-H'—V (Kawasaki) 
exchanges: one in the x and one in the y direction. This 
can be described by the generalized Kawasaki update 



(3) 



with probability p for attachment and probability q for 
detachment. 

We can also call the '+l'-s as particles and the '— l'-s 
as holes living on the base square lattice, thus an attach- 
ment/detachment update corresponds to a single step 
motion of an oriented dimer in the bisectrix direction 
of the x and y axes. We update the neighborhood of the 
sub-lattice points, which are the crossing-points of the 
dashed lines. In [ill [l2| we derived how this mapping 
connects the microscopic model to the KPZ equation and 
investigated the surface scaling numerically. Our best es- 
timates obtained by simulations up to sizes L = 2 15 for 
the two dimensional KPZ universality class is in agree- 
ment with the operator product expansion result [261 ] 



a 



0.395(5), ,3 = 0.245(5), z = 1.58(10) (4) 



within error margin (12J. 

Besides, presenting some more interesting results 
about the universal scaling behavior of KPZ in this paper 
we move further, and extend our mapping for describing 
more complex surface reactions. In particular, we inves- 
tigate models with surface diffusion processes, which are 
relevant in material science. We show the emergence of 
patterns and follow their coarsening dynamics. 



A. The simulations 

Although the bit-coded simulations are run on the un- 
derlying conserved lattice gas of size L x L, starting from 
hi.i = 1 we reconstruct the surface heights from the dif- 
ferences 



(T X (l, l)+y ]<Ty(i,k) 



1=1 



fe=l 



(5) 



at certain sampling times (t), selected with power-law 
increasing time steps and calculate its width: 



1 



L 



1/2 



(6) 



In the absence of any characteristic length, the surface is 
expected to follow Family- Vicsek scaling [13], when we 
start from a flat initial condition 

W(L,t) oc t , for t «t«t s (7) 
oc L a , for t » t s . (8) 
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Here, a is the roughness exponent for t >> t s when the 
correlation length has exceeded the linear system size L; 
and P is the surface growth exponent, which describes 
the time evolution for earlier (non-microscopic t » to) 
times. The dynamical exponent z can be expressed by 
the ratio 

z = a/p. (9) 

However in case of pattern formation multi-scaling is 
present in the system and the roughness exponent cal- 
culated in different window sizes is not constant and sat- 
isfies a different, anomalous scaling law (see [H, HI])- 

The morphology of pattern formation in experiments is 
usually followed by the measurement of some characteris- 
tic size, e.g. the wavelength in periodic structures, which 
evolves non-trivially with time. We can easily define and 
measure such quantity in our model. By following the 
up, or down slopes in the x direction we have strings of 
T'-s or ' — 1' of length Sk in the fc-th slice of the lattice 
gas in the y direction. We shall characterize patterns by 
calculating the (y) average of the longest Sk value 

L 

A = 1/L max(sfc) . (10) 

This characteristic length, corresponds to the longest x 
slope, or to the slowest x mode in the Fourier decom- 
position, will provide information about scaling of the 
wavelength in our analysis. 

By the simulations we apply periodic boundary con- 
ditions in both directions, and start from the flat space 
corresponding to a zig-zag configuration of the slopes (see 
FigHJ, therefore it has a small initial width W 2 (L,0) — 
1/4. Before the scaling analysis we always subtract this 
constant, being the leading-order correction, from the 
raw data. Averaging was done usually for 100 — 200 sam- 
ples for each parameter value. 

In practice each lattice site can be characterized by 
the 16 different local slope configurations, but we update 
it only when the condition ([3]) is satisfied. Furthermore 
due to the surface continuity not all configurations may 
occur and we can describe a lattice site by using only 2 
bits. This permits efficient storage management in the 
memory of the computer and large system sizes. The 
updates can be performed by logical operations either on 
multiple samples at once or on multiple (not overlapping) 
sites at once. Our bit-coded algorithm proved to be ~ 40 
times faster than the conventional FORTRAN 90 code. 
A crucial point is to use a good, high resolution random 
number generator, because in case of the p — 1 KPZ pro- 
cess the only source of randomness is the site selection, 
which must be done in a completely uniform way. Oth- 
erwise we realize a KPZ with quenched disorder which 
belongs to a different universality class (see 16]). Wc 
used the latest Mersenne- Twister generator [30] in gen- 
eral, which has very good statistical properties and which 
is very fast, especially by the SSE2 instructions. But we 



tested our results using other random number generators 
as well. 

An elementary Monte Carlo step (MCs) starts with a 
random site selection. This is followed by testing if the 
place is appropriate for update i.e. 'roof-top' for detach- 
ments or 'valley-bottom' for attachment ([3]). The update 
is done with the prescribed p and q probabilities and the 
time is incremented by 1/L 2 , such that one MCs corre- 
sponds to a full lattice update. Throughout the paper 
we use this unit of time. 



B. Generalizations of the octahedron model 

An obvious first step is to combine the deposition and 
the removal processes creating a conserved dynamics. 
A simultaneous octahedron detachment and deposition 
in the neighborhood can realize an elementary diffusion 
step . Surface diffusion is a much studied basic process 
j28{ . Several atomistic models have been constructed and 
investigated with the aim of realizing Mullins-Herring 
(MH) diffusion [3ll . [HJ and scaling (for a recent review 
see (33). 

The Langevin equation of MH is a linear one, with a 
V 4 lowest order gradient term 

d t h{x, t) = ^ 4 V 4 /i(x, t) + r/(x, t) . (11) 

emerging as the result of a curvature driven surface cur- 
rent j(x,t) cx V(V 2 /i(x, t)), which obeys the conserva- 
tion law 

&ft(x,t)+Vj(x,t) =r?(x,i) . (12) 

Here the noise rj (x, t) is a non-conserved, Gaussian white 
one, which can be the result of fluctuations in the ion- 
beam intensity directed against the surface. This equa- 
tion is exactly solvable, and exhibits a scaling invariance 
of the roughness characterized by the exponents: 

a=(4-d)/2, /3 = (4-d)/8, z = 4 (13) 

below the upper critical dimension d c = 4. 

Microscopic models realizing this behavior are mainly 
unrestricted solid on solid (SOS) type, that can pro- 
vide steep slopes and strong curvatures necessary for 
the a > 1 roughness exponent in one and two dimen- 
sions. However it turned out that the asymptotic uni- 
versality class of the various limited mobility growth 
models is a surprisingly subtle issue. Many of the ear- 
lier findings proved to be incorrect due to pathologi cally 
slow crossover and extremely long transient effects [34j . 
Anomalous scaling, by which the local and global be- 
havior is different, has been found to be relevant in such 
'sup er-rough' models, where large local slopes are present 
p9l [35l - [39| . In fact, according to our knowledge, only the 
"larger curvature (or Kim-Das Sarma) model" [40, |4jJ 
and the "n — 2 model" of [42j] exhibit MH universality 
class scaling asymptotically 



4 



The scaling in other 'atomistic' models crosses over to 
behavior dictated by more relevant terms in the sense 
of renormalization group. This can be the EW class if 
V 2 (/i) 43] is present, or the Molecular Beam Epitaxy 
(MBE) class in case of the fourth order non-linearity 
V 2 [V(/i) 2 ] 0-113] 0. The nonlinear MBE equation: 

a t /i = */ 4 V 4 /i + A2 2 V 2 [V(» 2 ]+r? (14) 

with non-conserved, Gaussian noise r\ is just the con- 
served version of KPZ (CKPZ) and exhibits the following 
scaling exponents: 

a=(4-d)/3, P = (4-d)/(8+d), z = (8+d)/3 . (15) 

As we can see all exponents are smaller than those of the 
linear MH (fT3|) class values and differ from those of the 
two-dimensional KPZ class ^ significantly. 

Here we present RSOS models with Ah = ±1 height re- 
striction within the framework of our previous approach 
[TT| , which in the limit of weak external noise exhibit MH 
or MBE scaling. Due to the simple construction these can 
be mapped onto lattice gases of diffusing dimers, allowing 
an easier way to study the effects of MH and MBE sub- 
processes of more complex system. Earlier less restrictive 
RSOS models were used to describe these classes espe- 
cially in one dimension (for a review see (33[). 

A further step in generalizing our surface models will 
be the combination of different sub-processes resulting in 
nonequilibrium system. For example, by adding a com- 
peting MH diffusion to the KPZ updates one can model 
the noisy Kuramoto-Sivashinsky (KS) equation [H, [49| . 

d t h = v + aV 2 h - v^h + \ 2 {Vh) 2 + r, . (16) 

Here the surface tension coefficient a is negative (in con- 
trast with the KPZ), whereas is a positive surface 
diffusion coefficient. However in our simulations we re- 
alize a gauge transformed situation: normal KPZ (with 
positive a) competing with an inverse MH (iMH) with 
negative (78|. In KS even the deterministic variant 
(j) = 0) exhibits spatio-temporal chaos, and is useful to 
describe pattern formation, such as chemical turbulence 
and flame- front propagation [48]. There are other phys- 
ical systems, including ion sputtering, where the noisy 
version of the KS equation is used [50] . In one dimension, 
field theory has proved that KS belongs to the 1 + 1 di- 
mensional KPZ universality class [50j. However in higher 
dimensions perturbative field theory cannot access the 
strong coupling fixed point. Numerical studies [Hl1 - F53| 
have provided controversial results in 2 + 1 dimensions, 
hence it remained a controversial and challenging prob- 
lem to clarify the asymptotic scaling behavior of KS [54l - 

H3- 

It was pointed out [5|| [59| that grooved phases and 
growth instabilities may emerge as the consequence of 
broken detailed balance condition: 

P({*})«W = P({i'}) Wi ,-x (17) 



where P({i'}) denotes the probability of the state {i} 
and Wi->i> is the transition rate between states {i} and 
{i 1 }. This means that complex structures and patterns 
can emerge in nonequilibrium system. In one dimension 
a model of massive particles exhibiting momentum has 
been shown to exhibit KS scaling behavior [6(j. Later 
another one-dimensional RSOS growth model was con- 
structed [6l| , in which deposition and diffusion of single 
ad-atoms were competing. It was suggested that large- 
scale behaviour could be described by the noisy KS equa- 
tion. 



II. REALIZING THE MBE SURFACE 
DIFFUSION 

As we have mentioned in the Introduction we gen- 
erate surface diffusion as the simultaneous adsorp- 
tion/desorption of octahedra. Therefore, in our model 
after an appropriate removal site (a roof-top) selection 
has been done, we search for a valley bottom place in the 
neighborhood for deposition. The target site is chosen in 
the ix or ±y direction, with the probabilities: p+ x , P-x 
or p+ y , P- y respectively (see FigllJ. Throughout of our 
studies we normalized the attempt probabilities. The 
maximal jump distance was fixed to be l m < 4 lattice 
units following computer experiments. According to the 
construction the nearest neighbour jump, corresponding 
to intralayer diffusion, requires l m = 2, while larger jump 
sizes allow intra-layer transport. To create MBE or MH 
behavior we must allow inter-layer transport. Consid- 
ering jumps with l m > 3 does not make difference in 
the scaling, larger jumps describe faster diffusion on the 
expense of more CPU time. Having no apriori knowl- 
edge of the distance dependence of jump sizes and ex- 
pecting insensitivity of universal scaling on the rates of 
short ranged interactions we used a constant probability 
for each direction (for more discussion see Sect. Ill A"| . 

To control this kind of surface diffusion we impose ad- 
ditional constraints for the accepting a move. We have 
tried two kinds of rules based on the local neighborhood 
configurations. The first one is very simple and requires 
that the height of a particle at the final state, is higher 
than that of its initial site 

hfi„ - h ini > , (18) 

which makes a surface more rough (inverse/roughening 
diffusion) . The second one is based on the local curvature 
at the update sites as will be discussed in Sect. Ill Bl 

A. The larger height octahedron diffusion model 

By this simple condition, we expect to generate rough 
surfaces of large curvatures, since an octahedron jumping 
to a higher position has usually smaller number of neigh- 
bours in the lateral direction (in this model we do not al- 
low particles to evaporate from the bottom of a V shaped 
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valley). In the language of lattice gases regions of large 
curvature and maximal slope correspond to large regions 
of dense particles and holes. An attractive interaction of 
the dimer moves enhances such segregation. Contrary, 
we can realize a smoothing/normal diffusion, when we 
accept jumps with the condition: hfi n — hi n i < 0. 

By increasing the local heights the formation of 
pyramid-like structures , unstable growth, occurs simi- 
larly as in the n = 2 SOS model [42j, but in our case 
the slopes are limited to 45 degree. However, by appro- 
priate length rescaling any sharp, continuous surface can 
be approximated. In the forthcoming, we will investigate 
both the inverse and the normal diffusion cases and call 
it the larger height octahedron diffusion (LHOD) model. 

Starting from a zig-zag initial condition, corresponding 
to the flat surface, the hopping with the condition (fl"8|) . 
one can make the surface more rough, however after the 
slopes of size l m are developed the evolution stops, be- 
cause the octahedra are not allowed to pass longer 45 
degree gradient sections (see Fig. [5]). To overcome this, 
we have been trying to allow arbitrarily long jump sizes 
(non-local model), with different, heuristic jump distance 
dependent probability functions. We did not find an ap- 



Initial state Frozen state 




• O • O • G • G - • G G • • • G G 



FIG. 2: (Color online) One dimensional view of a diffusion 
hop of octahedra by 3 lattice units to the right. The slopes 
mapped onto the lattice gas shown below the surface. This 
roughening surface diffusion move corresponds to two simul- 
taneous, attracting Kawasaki exchanges of the gas (dimers in 
two dimensions). Starting from the flat (zig-zag) initial state 
of the pure octahedron model the (|18|l process freezes follow- 
ing slopes of maximal length are developed, which can't be 
over-jumped. 



propriate one that could produce the expected MH scal- 
ing behavior, instead we realized Levy flight like mod- 
els with anomalous diffusion 62] exhibiting non-universal 
scaling. In he light of recent field theoretical interest [63] 
these can also be the target of further investigations. 

However, non-local models are rather complex and con- 
nection to reality is not always straightforward. There- 
fore, we followed an other strategy by adding a small 
amount of extra randomness of EW type to our short- 
range, binary RSOS model. This means the addition 
of random adsorption/removal events with small prob- 
ability p = q « 1 among the LHOD updates. Such 
events can break up the barriers by splitting up the long 
monotonous slopes built by the LHOD dynamics. If they 
are done very rarely, i.e. with less than 100 times smaller 



probability than the diffusion attempts, they can influ- 
ence the very late asymptotic scaling behavior only, caus- 
ing an ultimate crossover to the EW class (see for exam- 
ple [sl). 

In reality, and in the models we are about to study, 
such randomizing effects are always present, thus we are 
satisfied, if we can confirm numerically the MBE scal- 
ing for intermediate times. It was not our principal aim 
to provide a model, which exhibits pure MBE type of 
scaling in the thermodynamic limit. Therefore we have 
performed computer experiments on these models as dis- 
cussed below. 




icr 5 icr 4 icr 3 icr z icr' 10° io~ 5 icr" icr 3 icr 2 icr' 10° 
t/L 4 t/L 4 

FIG. 3: (Color online) Data collapse of the anisotropic LHOD 
model assuming MH class exponents for p+ x — 1 (diffusion 
to the right) with small EW noise p — q — 0.01 (a), p = q = 
0.005 (b), for sizes L = 32,64, 128 (top to bottom at the left 
side). For the growth exponent fitting (dashed line) results in 
/3 = 0.26(1). 



We found that the addition of the small (EW) noise 
(p = q « 1) sustains the surface currents, and in case 
of spatially anisotropic diffusion (p+ x = 1, p- x — P y — 
P- y = 0) the scaling becomes Mullins-Herring type (see 
Fig. EJa)), characterized by the exponents a = 1, /3 = 
1/4, z = 4 in two dimensions. The collapse of curves 
is very good in the vertical direction, corresponding to 
a = 1 , and the horizontal scaling improves as we decrease 
the noise (see Fig. [3Jb)). In opposite, by increasing the 
aplitude of the EW noise, we can find better collapse with 
smaller dynamical exponent, describing the cross-over to 
the EW behavior, which has z = 2. The time dependence 
shows deviations from the pure scaling-law, especially for 
early times, still before the saturation the growth of W(t) 
can be fitted with the exponent f3 = 0.26(1). [79j This 
suggests that in the p = q — > limit the true MH scaling 
emerges. Simulating larger sizes is very hard, because 
due to the large dynamical exponent z the saturation is 
shifted to very late times (for L = 128 this happens for 
t > 2 x 10 8 MCs only). 

In case of isotropic diffusion the scaling is MBE class 
type in general (see Fig. 0|a) for l m =4), characterized 
by the exponents a — 2/3, (3 — 0.2, z — 10/3 in two di- 
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mensions. Therefore, the algorithm with the LHOD up- 
date (|18p breaks the detailed balance condition (JTTJ) and 
introduces a non-linearity. However, this non-linearity is 
small and by decreasing l m it becomes even smaller (see 
Fig. |4] (b)). One can find a rather good collapse with 
the MH exponents in case of l m — 3. We estimated the 
growth exponent by calculating the local slopes 



&//(*) = 



In W(t, L oo) - In W(t', L ->• oo) 
ln(i) - ln(i') 



(19) 



As one can read-off from the inserts of Fig. [U for l m = 4 
this effective exponent extrapolates to the MBE value 
{fi = 0.20(2)), while for l m = 3 to = 0.26(1) agree- 
ing with the MH exponent in the t — ¥ oo limit. This 




of squares of the projected octahedra. As Figure [5] shows 
one can describe the curvatures c x (i,j) (x £ (x,y)) by 
the products (or differences) of the local slopes 



(20) 



At each update we calculate the sum of the change of 
local curvatures at the origin and the target 
sites 

AH = A E E c x(^')+ A E E c *(^")< ( 21 ) 

X=x,y(i,j) X=x,V (i'./j') 

where (} denotes the plaquette neighborhood sites as 
shown on Fig. [5j This gives maximal value: H — 4 for a 
local tops and bottoms and the minimal value: H = —4 
for a locally flat (zig-zag) configuration. Using this value 

Y x 

\0 y (i+i,j+2) a x (i+2,jy y 

c >( i+1 <i +1 t ■ C£+2j) 




<?x(i-lj) 



<^i+lj-2) 



FIG. 4: (Color online) Data collapse of the isotropic LHOD 
model assuming MBE (left, l m — 4) and MH (right, l m — 
3) class exponents in the presence of small EW noise p — 
q = 0.005 in sizes L — 32, 64, 128 (top to bottom at the 
left side). For the growth exponent power-law fitting (dashed 
line) results in /3 = 0.26(1). The inserts show the effective ft 
exponents for L = 64. 

means that the LHOD rule introduces more possibilities 
for breaking the detailed balance condition (|17p . when 
more degrees of freedom (more directions or larger jumps 
sizes) are allowed. This provides an explanation for the 
difference between the scaling behavior of the isotropic 
and anisotropic diffusion. To avoid such non-linearity 
completely we introduced a more complex update rule, 
based on the local curvature conditions keeping the spa- 
tial symmetries. 

B. The larger curvature octahedron diffusion 
model 

In this section we describe a transition probability in 
addition to the octahedron surface hopping model, which 
satisfies the detailed balance condition (TTTf , thus enables 
one to realize linear, equilibrium MH diffusion steps. The 
local curvature of the surface is calculated at the 4 edges 



FIG. 5: Local slope a x (i,j) (circles) and curvature c x (i,j) 
(squares) variables at an update site. Filled circles correspond 
to upward, empty ones to downward slopes of the surface. 
This plaquette configuration models a valley bottom site, with 
the total curvature H — 4. 

we accept the update with an Arrhenius type of proba- 
bility 

Wi-n< = 1/2 [1 - atanh(-A J ff 2 )] , (22) 

where a is a constant. This form is very similar to what 
was used in case of the one-dimensional n — 2 model 
[42| and enhances (suppresses) roughening moves if a > 
(a < 0) respectively. In [42J symmetry arguments was ap- 
plied for the lowest order series expansion of of the 
model to prove a connection with the MH equation. Now 
similar derivation can be done by extending the model for 
dimer variables in two dimensions. In the forthcoming we 
shall call this the larger curvature octahedron diffusion 
(LCOD) model. 

We have simulated the LCOD with the parameters: 
a = 0.1, l m = 3, corresponding to inverse/roughening 
surface diffusion and p — q — 0.05. The scaling behavior 
has been found to agree very well with that of MH univer- 
sality class values even in case of isotropic diffusion (see 
Figure[5]). The effective j3 e ff(t) converges to f3 — 0.25(1) 
before the saturation. The wavelength grows logarithmi- 
cally in time (see insert of Figure [5]) and after a steady 
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A. Spatially anisotropic surface diffusion 

From the point of pattern formation the LHOD and 
LCOD model behaves differently. We always start the 
simulations from a flat surface and watch if stable pat- 
terns can arise. In case of an anisotropic inverse LHOD 
model of diffusion probability p± x = 1, p± y = a com- 
peting EW process always generates ripple patterns as 
shown on Fig. [7]), which is stable for all p = q < I. This 
formation is metastable against KPZ (height anisotropy), 
but for very large times (in the steady state) the rip- 
ples become uneven, blurred and cut into smaller pieces. 
The wavelength, defined as (| 10[) grows only a little (in a 



FIG. 6: (Color online) Scaling behavior of the isotropic LCOD 
model for L = 32, 64, 128 (top to bottom at the left side). The 
data collapse has been achieved with the MH class exponents. 
For the growth exponent fitting (dashed line) results in ft = 
0.25(f). The insert on the right shows the same by local 
slopes. The insert on the left shows the evolution of A. 



state value is reached it scales logarithmically with the 
the system size too. 

For anisotropic diffusion (p± x — l,p±y — 0) we have 
obtained similarly good MH class surface scaling, with 
logarithmic time and size dependence of A again. 

For a = 0, without any EW noise, one can find loga- 
rithmic growth in the LCOD model in time 

W(t, L ->■ oo) ocln(t) (23) 

and logarithmic surface roughness dependence of 



W(t -> oo, L) oc ln(L) 



(24) 



Data collapse fitting for the dynamical exponent on the 
other hand results in space-time anisotropy with z = 
4. This means that for the a — 0, noiseless case we 
could realize the universality class behavior of MH with 
conserved (purely diffusive) noise, characterized by the 
exponents: a = /3 = 0, z = 4 (see |16|). 



III. PATTERN GENERATION BY COMPETING 
INVERSE MH AND KPZ PROCESSES 

As we have shown in the previous section, in the zero 
noise limit our RSOS surface diffusion processes gener- 
ate growth with MBE or MH scaling behavior. In this 
section we investigate them in the presence of compet- 
ing KPZ updates. In our simulations we initiate hopping 
with probabilities: p+ x , P-x , P+y, P-y alternately with 
the deposition (with probability p) and removal (with 
probability q) processes. We follow the surface rough- 
ness and pattern formation with the corresponding wave- 
length growth. 




100 200 300 400 



FIG. 7: (Color online) Snapshot of surface heights of the 
ripple patterns generated by the parameters p± x — 1, 
p± y = (anisotropic, inverse MH) and p = q = 1 adsorp- 
tion/desorption at t = I0 4 MCs, in the LHOD model of linear 
size L = 512. 



power-law manner) and saturates quickly. 

However, if we create such anisotropy in which a 
steady, direct current (DC) flows, for example when only 
p +x = 1 and all the others are zero we can find a different 
behavior. In this case, for weak EW or KPZ the ripples 
are not completely straight and exhibit a coarsening as 



A oc i a24 « 



(25) 



The ripples are more straight in the KPZ case than for 
up-down isotropic deposition/removal. Furthermore, a 
good data collapse can be obtained with the MH class 
exponents for sizes L = 32, 64, 128 as shown on Fig. [8] 
However, this scaling can be destroyed by increasing the 
strength of the non-conserved reaction and the power- 
law crosses over to logarithmic growth of A. fn case of 
strong KPZ deposition (p = 1, q = 0) the asymptotic 
scaling of the LHOD becomes completely KPZ type, the 
wavelength saturates very quickly (see Fig. [S]) and the 
ripple structure smears. In principle we should expect 
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FIG. 8: (Color online) The wavelength growth in the LHOD 
model for anisotropic diffusion with steady DC current p+ y — 
1, p = q = 0.005 for sizes L = 32, 64, 128 (top to bottom at 
the beginning). Dashed line: power-law fit with the exponent 
f3 — 0.24(1). The left insert shows the corresponding pattern. 
The right insert corresponds the isotropic diffusion case p± x — 
p± y — 1, where X(t) grows logarithmically. 
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FIG. 9: (Color online) Data collapse for deposition (p = 
1, q = 0) and anisotropic, inverse diffusion p± y — 1 
in the LHOD model with KPZ class exponents for L — 
64, 128, 256, 512, 1024 (top to bottom curves at the right side). 
Right insert: A(r) for L — 1024, left insert the blurred ripple 
structure. 



anisotropic KPZ behavior here, but it is well known that 
in two dimensions such spatial anisotropy is irrelevant, 
hence the isotropic KPZ class behavior |65j is not sur- 
prising. Note that our anisotropic KS model is different 
from what is called and considered to be the " anisotropic 
KS" in the literature [66[ , because in our case the surface 
diffusion (corresponding to the V 4 term) is anisotropic. 
Such models are very hard for analytic treatment and 
spatial anisotropy is introduced in the V 2 terms usually. 
The anisotropic LCOD is less effective for ripple for- 



mation than the LHOD. In this case the patterns are 
smoother, and the A scales logarithmically in time and 
by the size. The only exception is when we allow a steady 
DC current again. In this case similar power laws as 
shown on Fig. [9] emerge for weak EW or KPZ. 

When DC current is not allowed (p± x = l,p± v = 0), 
only spatial anisotropy in the LCOD and we add a KPZ 
(p = 1, q = 0) update we can see the emergence of KPZ 
scaling. In this case the wavelength depends logarith- 
mically both on time and the size L. It is important 
to realize that for short, one decade length, time win- 
dows the wavelength growth can also be well fitted with 
a power-law 



X(t,L 



(26) 



which resembles to experimental results, but since the 
steady state values exhibit a clear logarithmic depen- 
dence on the sizes 



, L) oc ln(L) 



(27) 



we don't think that this 'power-law' fit would corre- 
spond to a real asymptotic behavior in the thermody- 
namic limit. 



B. Spatially isotropic surface diffusion 

When the isotropic, inverse surface diffusion competes 
with the (smoothing) EW process one can observe dot 
formation both in the LHOD and LCOD models. Fig.fTUl 
shows a snapshot of the growing dots for LHOD in the 
presence of weak EW process. Here, one can see rectan- 
gular shaped patterns corresponding to the lattice sym- 
metry. In case of LCOD the contrast of the patterns is 
smoother, roughly circular. However, this pattern coars- 
ening is much slower than in case of ripples. 

We shall discuss the LHOD results first, which accord- 
ing to our previous numerical analysis, corresponds to 
the nonlinear equation (KPZ+MBE): 

d t h = aV 2 h + X 2 (yh) 2 + iy 4 \7 4 h + X 2 2V 2 [V(h) 2 ]+ri (28) 

In case of EW type of deposition/removal the nonlinear 
term vanishes (A2 = 0) and in fact we model the CKPZ 
behavior. The characteristic size (A) of the dots grows 
logarithmically (see insert of Fig. [5J in time and X max oc 
ln(L). The pattern formation is more pronounced in the 
LHOD case than in the LCOD model. We associate it to 
the up/down anisotropy present in the LHOD. 

The same kind of patterns can be also observed in case 
of strong KPZ anisotropy p = 1, q = for short times, 
but later the dots are smeared. For the LHOD+KPZ case 
a very slow crossover to KPZ scaling (see Fig.QT]) occurs. 
Although the data collapse with KPZ exponents is rather 
poor for smaller sizes (it would be better with larger z 
and a exponents corresponding to the MBE class), for 
L > 512 it agrees with KPZ. The wavelength saturates 
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precision KPZ simulation result [12]. The wavelength 
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FIG. 10: (Color online) Snapshot of surface heights of the 
dot patterns generated by the LHOD model with parameters: 
P±x = p± y = 1 (isotropic, inverse MH) and p = q = 0.1 at 
t = 10 4 MCs. 



FIG. 12: (Color online) Data collapse of the L = 128, ...1024 
(top to bottom at the right) LCOD model (p± x = p± y = 1) 
with competing deposition (p = 1). One can see clear KPZ 
scaling. The left insert shows the logarithmic growth of A for 
L = 128,256,512 (bottom to top). Right insert: /3 e // as the 
function of time. 




FIG. 11: (Color online) Data collapse of the L = 128, ...1024 
LHOD model (p± x = p± y = 1) with a competing deposi- 
tion (p = 1) process. One can see a very slow crossover to- 
wards KPZ scaling. The right insert shows the growth of A 
for L — 512. The left insert is a snapshot of the steady state, 
corresponding to the smeared KPZ height distribution. 



very quickly (see insert of Fig. Illj) in agreement with 
KPZ, where no coarsening is expected. 

Having confirmed that the LCOD model exhibits MH 
scaling, now we can investigate the scaling behavior of 
the (inverse) KS equation fp~6|) . described by the combi- 
nation of inverse MH and normal KPZ processes. We 
have run extensive simulations up to t = 3 x 10 6 MCs 
(for L = 128, 256, 512, 1024) to obtain firm numerical ev- 
idence. As Fig. [T2] shows the finite size scaling collapse 
with with KPZ exponents is satisfied and the effective /3 
extrapolates to 1/4. This value agrees well with our high 



grows logarithmically in time (|26[) (see insert of Fig. [T2]) 
and saturates well before the steady state. In the steady 
state it grows slowly with the system size as (|27p. In a 
one decade long time window one can fit the data with 
X(t,L — > oo ) cx £°- 12 ( 2 ) ; but due to the clear logarithmic 
behavior in the steady state (|27p. one should not take 
such power-law fitting very seriously. 

For the sake of completeness we have performed sim- 
ilar analysis for the anisotropic LCOD model (p± x = 
l,p± y = 0) with KPZ too, for sizes L < 1024. We have 
obtained agreement with the KPZ scaling for the width 
W and logarithmic growth for A as in case of isotropy. 



IV. KPZ IN THE PRESENCE OF NORMAL 
SURFACE DIFFUSION 

In this part we show the scaling behavior study of 
the KPZ process in the presence of normal (smooth- 
ing) LHOD, introduced in the previous section. First 
let's consider the weak isotropic diffusion case: p± x = 
P±y = 0.1. As Fig. [TBI shows the W(t) curves for sizes 
L = 64, 128, ...2048 exhibit a good collapse when we 
rescale them with the 2+1 dimensional KPZ exponents 
(J4j) . This means that the KPZ scaling is stable against 
of the introduction a smoothing, MBE type of surface 
diffusion. 

On the other hand when we add strong LHOD diffu- 
sion, p± x = p± y = 0.9, to the KPZ process (p = 1) the 
surface growth slows down and we don't find the KPZ 
scaling any more (see the lower part of Fig. [T3]) . In- 
stead, logarithmic surface growth emerges, as shown in 
the insert of Fig. [T3J A fitting with the form W 2 (t) = 
A + B ln(i) for the asymptotic growth regime gives A — 
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FIG. 13: (Color online) Data collapse of KPZ deposition (p = 
1) and weak, isotropic normal LHOD (higher curves) for L — 
64, 128, ...2048 (top to bottom). In case of strong diffusion 
(lower curves) the KPZ scaling disappears and as the insert 
shows logarithmic growth can be observed. 



0.40(1) and B = 0.26(1). The amplitude of this growth 
is different from the exactly known universal value of the 
EW class in two dimensions: A E w — 0.151981 [67j . 

The wavelength saturates very quickly to the maxi- 
mal value, which for weak diffusion scales logarithmically 
with the system size (|27p. while in case of the strong 
diffusion the characteristic length remains on the order 
of lattice unit: \ ma x — 1, with a very weak size de- 
pendence, corresponding to uncorrelated surface heights 
(see Fig. IT4"|) . It is well known, that in the strong dif- 
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FIG. 14: (Color online) The wavelength saturates quickly for 
KPZ + weak LHOD (higher curve) and KPZ + strong LHOD 
(lower curve) diffusion (L — 2048). The insert shows \ ma x 
versus L. 



fusion limit, the relevant fluctuations below d c can be 
washed away, hence the logarithmic surface growth we 
observe should correspond to the mean-field behavior of 
the strong-coupled KPZ. This provides us a unique way 



to study the crossover behavior between the KPZ and 
KPZ mean-field behavior. 



V. PROBABILITY DISTRIBUTION RESULTS 

Here we present the probability distributions P(W 2 ) 
obtained in the steady state of our model. Such distri- 
butions are universal, hence they complement the previ- 
ous scaling results. The exact functional form for KPZ 
is known in one dimension only, but in two dimensions 
very precise numerical data exist, obtained via other sur- 
face models [68|. The distributions of those KPZ models 
have been determined in higher dimensions [6!| , suggest- 
ing the lack of finite upper critical dimension. 

First, we compare our P(W 2 ) results for KPZ with 
those of [iil in d — 2, 3, 4, 5 dimensions. The W 2 distri- 
bution data were taken from the saturation regimes and 
analyzed in systems of sizes: L — 1024 (2d), L = 512 
(3d), L = 64 (4d), L = 32 (5d). The presented data are 
coming from our higher dimensional KPZ simulations, 
using the extended octahedron model described in [il^ . 
When we rescaled our data with (W 2 ) as Fig. [Pol shows 
we found very good agreement in d = 2,3,4,5 dimen- 
sions with the earlier KPZ distribution curves. Again, 




w 2 /<w 2 > 



FIG. 15: (Color online) Comparison of the P(W 2 ) of the 
higher dimensional octahedron model results (symbols) with 
those of 69] (lines) in d = 2, 3, 4, 5 spatial dimensions (bottom 
to top). 



we can't see a signal for an upper crictical dimension at 
d c = 4, conjectued by theoretical approaches (see for ex- 
ample (70j|). 

Furthermore, we tested our surface scaling results for 
KPZ with in the presence of diffusion. In case of com- 
peting KPZ and LHOD/LCOD processes (i.e. for p = 1, 
9 = 0, p± x = p± y = 1) we determined the P(W 2 ) dis- 
tributions well in the saturation regime of systems of 
size L = 1024. The steady state could be reached for: 
t >~ 10 5 MCs in case of KPZ+LHOD and for t > 5 x 10 7 
MCs in case of the KPZ+LCOD model. We generated 
100 independent samples, and cut out the steady state 
data of W 2 (t). We calculated the scaled steady state 
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probability distributions as shown on Fig. [16] for differ- 
ent combinations. The comparison with the histogram 
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FIG. 16: (Color online) Comparison of P(W 2 ) of the 
KPZ+LHOD (black boxes); KPZ+inverse LHOD (blue 
dots); KPZ+inverse, anisotropic LHOD (pink rhombuses); 
KPZ+inverse LCOD (orange triangles) with that of the KPZ 
from ref. [69l] (solid line). 

of the 2 + 1 dimensional KPZ from 69] shows a good 
agreement in general, providing further numerical ev- 
idence that the KS model asymptotically exhibits the 
KPZ universality class scaling. Our results complement 
the one-dimensional simulations results of [611 ] . in which 
the equivalence of KS and KPZ scaling was confirmed 
numerically. 

VI. CONCLUSIONS AND OUTLOOK 

We have returned to some unresolved questions of ba- 
sic surface growth phenomena using extensive computer 
simulations of atomistic models. Contrary to earlier stud- 
ies we could perform numerical analysis of models with 
surface diffusion in 2 + 1 dimensions, due to the effective 
mapping of RSOS models onto binary lattice gases. 

We have shown that in the zero external noise limit 
RSOS models with short range interactions can be con- 
structed, which exhibit Molecular Beam Epitaxy or 
Mullins Herring type of surface growth. For inverse 
(roughening) diffusion, which increases the local curva- 
ture, unstable growth, resulting in pyramid-like struc- 
tures emerges. The size of these structures is limited 
only by L, which is not directly comparable with real ma- 
terials. We created these MBE or MH type of atomistic 
models in order to study them in a competition with non- 
conserved KPZ processes. The simulations provided nu- 
merical evidence that strong, smoothing surface diffusion 
can slow down the KPZ to a logarithmic growth, thus we 
are able to reach the mean-field behavior of the strongly 
coupled KPZ fixed point in two dimensions, which is ex- 
pected to show up in high dimensions only. 

The mapping of the surface models onto lattice gases 
implies that the (anisotropic), oriented diffusion of 



dimers (KPZ) is stable against the introduction of an 
attracting force among them, but a strong repulsion can 
destroy the fluctuations, resulting in a mean-field behav- 
ior. We provided strong numerical evidence using surface 
scaling and probability distribution studies, that the KS 
model exhibits KPZ scaling in 2 + 1 dimensions as con- 
jectured by field theory. We summarized the models we 
considered in case of spatially isotropic surface diffusion 
and non-conserved noise in Table |U 

Further studies of LCOD, with different boundary con- 
ditions can set the target of the research of the surface 
inclination by mapping the surface tilt onto the total par- 
ticle concentration of the lattice gas. In particular the 
angle dependence of the phase transitions among differ- 
ent growth phases can be understood by considering the 
underlying driven gas. Using our method one can trans- 
form results of disorder or of anomalous diffusion between 
the surface and lattice gas models. 

We introduced a characteristic length scale A to fol- 
low the dynamics of patterns, which occur, if normal 
(smoothing) KPZ competes with inverse (roughening) 
diffusion. We investigated this race for MH and MBE 
process, with and without spatial anisotropies. In case 
of uni-axial surface diffusion ripple, while for x/y lattice 
isotropy, dot like pattern formation could be achieved. 
The wavelength growth is slow and saturates much be- 
fore the steady state. Usually, we found logarithmic time 
and system size dependence of A, except when steady DC 
current flows through the system. In this case the inter- 
faces are more rough, the ripples are bended and power- 
law scaling is observable. In this case the scaling of A 
agrees with the scaling of the width, i.e. characterized 
by the MH class exponents. This finding agrees with 3d 
kinetic Monte Carlo simulations of epitaxial growth and 
erosion on (110) crystal surfaces [7l[ and with analytic 
arguments [721 ] . Furthermore in case of ion beams with 
grazing incidence the dislocation dynamics results in such 
athermal, kinetic coarsening of the patterns [73| . 

The wavelength behavior can be understood with the 
help of considering the underlying lattice gas model, since 
the ripple or dot structures correspond to a phase sepa- 
ration. The coalescing surface pattern dynamics can be 
mapped onto the generalization of the reaction-diffusion 
process [l6| of extended objects. It was shown that in 
ASEP type of models, where strong phase separation is 
present the domain growth follows slow, logarithmic be- 
havior [73, IzU in case of smooth surfaces. On the other 
hand for rough surfaces, power-law coarsening of A has 
been derived using simple scaling arguments [76] . 

Finally we mention that our models enable effective, 
bit-coded, stochastic cellular automaton type of simula- 
tion of surfaces, hence they could be run extremely fast 
on advanced, graphic cards. 
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TABLE I: Overview of the surface models with spatial isotropy and non-conserved noise. Columns 2-5 show the sign of couplings 
of the differential equations considered. 
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